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1.0  SUMMARY 


This  expository  report  discusses  fundamental  aspects  of  the  polynomial  chaos  method  for 
representing  the  properties  of  second  order  stochastic  processes.  As  originally  developed  by 
Norbert  Weiner,  a  polynomial  chaos  represents  key  properties  of  a  stochastic  process  through  the 
application  of  finite  series  of  orthogonal  polynomials.  The  attendant  polynomial  expansion  is 
used  to  describe  the  statistical  properties  of  a  stochastic  process  based  upon  an  input  uncertainty. 
The  statistics  of  a  random  process  is  given  by  evaluating  the  appropriate  polynomial  chaos  for  an 
input  uncertainty  represented  by  one  or  more  random  variables.  An  evolved  application  of  this 
idea  applies  a  polynomial  chaos  to  represent  uncertainties  in  boundary  or  initial  conditions  for 
partial  differential  equations.  Here,  the  elementary  theory  of  the  polynomial  chaos  is  presented 
followed  by  the  details  of  a  number  of  example  calculations  where  the  statistical  mean  and 
standard  deviation  are  compared  against  exact  solutions.  The  Legendre  chaos  is  described  in 
some  detail  for  unifonnly  distributed  input  random  variables.  Also,  the  Hermite  chaos  is 
discussed  for  random  variables  possessing  a  Gaussian  distribution.  A  number  of  example 
problems  are  solved. 
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2.0  INTRODUCTION 


Real  world  problems  involving  science  and  engineering  always  involve  uncertainty. 
Specific  uncertainties  in  the  form  of  imprecisely  known  input  parameters  propagate  through  a 
physical  mechanism  and  create  uncertainty  in  the  system  output  forming  a  stochastic  process 
involving  one  or  more  random  variables.  A  simple  example  considered  later  in  this  report 
addresses  the  oblique  shock  wave  formed  by  the  supersonic  flow  over  a  wedge.  In  this  scenario, 
we  presume  that  the  wedge  angle,  with  respect  to  the  freestream  direction,  is  imprecisely  known. 
The  wedge  angle  is  regarded  as  a  random  variable  with  a  mean  value  and  a  standard  deviation 
based  upon  a  distribution.  The  ensuing  oblique  shock  wave  angle  becomes  a  function  of  this 
random  variable  and  takes  on  a  randomness  that  can  be  assessed  by  the  polynomial  chaos.  Other 
properties  associated  with  the  shocked  flow  field  also  adopt  random  properties  that  may  also  be 
assessed  by  the  chaos.  This  example  is  conceptually  simple  to  understand,  but  the  same  idea  can 
be  extended  to  other  more  complicated  physical  systems.  The  solutions  of  either  ordinary  or 
partial  differential  equations  subject  to  uncertainties  in  initial  or  boundary  conditions  also  inherit 
random  characteristics  that  can  be  revealed  by  the  chaos.  For  this  reason,  there  are  many 
applications  for  this  mathematical  tool.  Before  beginning  a  detailed  mathematical  development 
of  the  chaos,  it  is  instructive  to  discuss  the  nature  of  physical  uncertainties. 

2.1  The  Concept  of  Uncertainty 

As  is  the  case  for  physical  properties,  the  concept  of  uncertainty  has  been  analyzed  and 
properly  defined.[l]  Consider  the  uncertainty  existing  in  the  measurement  or  calculation  of  a 
physical  quantity.  If  the  quantity  possesses  a  true  value  with  no  variation,  then  the  variation 
sensed  in  measuring  or  calculating  the  quantity  is  referred  to  epistemic.  On  the  other  hand,  if  a 
physical  quantity  exhibits  a  natural  variability,  (that  is,  it  possesses  no  true  (or  certain)  value), 
then  the  associated  uncertainty  is  denoted  as  aleatoric.  In  the  Bayesian  probabilistic  framework, 
epistemic  uncertainties  can  be  addressed  by  techniques  such  as  those  described  below.  Aleatoric 
uncertainties  are  readily  treated  by  probabilistic  methods.  The  primary  application  considered  in 
the  results  section  of  this  report  can  be  regarded  as  having  epistemic  uncertainty.  In  this  case,  the 
wedge  angle  is  assumed  to  possess  a  small  variation.  With  the  use  of  clever  mathematical 
techniques,  the  statistical  parameters  (e.g.,  mean  and  standard  deviation)  of  functions  of  these 
uncertainties  can  be  predicted  with  accuracy. 

2.2  Mathematical  Representation  of  Stochastic  Processes 

The  mathematical  development  of  stochastic  processes  is  extensive,  the  culmination  of 
decades  of  research,  so  only  the  briefest  description  is  provided  here.  The  first  requisite  concept 
is  that  of  the  sample  space.  A  sample  space  is  the  set  of  all  possible  outcomes  of  an  experiment 
whether  the  experiment  is  simple  such  as  a  coin  toss  or  complex,  or  the  measurement  of 
temperature  or  fluid  stress  taken  at  some  point  in  a  turbulent  flow  field.  Detailed  examples  of 
sample  spaces  are  contained  in  Ross. [2]  On  a  given  sample  space,  much  like  a  mathematical 
domain,  we  define  functions  denoted  as  random  variables.  For  example,  on  the  event  space  for  a 
coin  toss,  we  can  define  a  random  variable  that  counts  the  number  of  heads  occurring  for  a 
certain  number  of  coin  tosses.  As  another  example,  we  may  define  a  random  variable  that 
represents  the  time  varying  electrical  voltage  existing  at  a  point  in  an  electronic  circuit.  The 
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former  is  a  discrete  random  variable  while  the  latter  is  a  continuous  random  variable.  A  discrete 
random  variable  takes  on  either  a  finite  or  countably  infinite  number  of  values  while  the 
continuous  random  variable  takes  on  a  finite  number  of  values.  Building  upon  these  concepts,  a 
stochastic  process  is  a  collection  of  random  variables  X{t)  indexed  by  a  set  T  where  t  e  T.  T  is 
denoted  as  the  indexing  set;  it  may  consist  of  either  a  finite  or  infinite  number  of  elements 
depending  on  the  nature  of  the  stochastic  process.  More  specifically,  X{t)  is  the  state  of  the 
stochastic  process  X  at  time  t.[ 2]  In  the  context  of  this  report,  X(t)  may  be  regarded  as  the 
uncertainty  in  some  quantity.  With  this  assertion,  we  may  envision  functions  of  the  uncertainty 
represented  by  X{t).  This  function  of  the  uncertainty  (random  variable)  requires  the  use  of 
“functionals”;  a  functional  may  be  thought  of  as  a  function  of  functions. [3]  In  this  sense,  a 
functional  is  defined  on  a  vector  space  of  functions,  and  this  concept  represents  a  significant 
mathematical  abstraction  beyond  that  of  the  more  commonly  known  vectors  defined  in  Euclidean 
space  9T . 

A  significant  difficulty  in  simulating  stochastic  processes,  or  more  particularly,  in 
simulating  stochastic  variation  in  otherwise  deterministic  systems  lies  in  the  lack  of  an  intuitive 
understanding  of  function  spaces. [4]  For  stochastic  problems,  the  function  space  must  be 
“measurable”  with  respect  to  a  probability  space  consisting  of  the  sample  (or  event)  space,  a  o- 
algebra  defined  on  the  event  space  and  a  probability  measure  P.  The  a-algebra  is  a  collection  of 
subsets  of  the  event  space  that  possess  special  properties. [5]  One  approach  to  the  problem  of 
simulating  random  or  stochastic  processes  is  Monte  Carlo  simulation,  a  technique  that  requires 
the  repeated  sampling  of  the  a-algebra.  This  method,  although  highly  effective,  requires 
sampling  many,  many  points  requiring  a  great  deal  of  computing  resources. [6, 7, 8]  Another 
approach  that  can  be  less  resource  intensive  entails  representing  the  stochastic  process  in  tenns 
of  a  set  of  orthogonal  functions  defined  on  the  function  space. [4,5]  Although  the  resulting 
orthogonal  decomposition  (a  spectral  approach)  consists  of  an  infinite  series,  the  series  can  be 
finitely  truncated  for  computational  purposes.  The  principal  function  space  containing  the 
orthogonal  function  set  is  denoted  as  L2(Q,P) where  Q  is  the  sample  space,  and  Pis  the 
probability  measure.  L2  is  a  Hilbert  space  known  as  the  space  of  square  integrable  functions 
widely  applied  throughout  mathematical  physics. [9]  For  a  given  stochastic  process  defined  at 
points  in  space,  an  orthogonal  decomposition  of  great  theoretical  interest  is  the  Karhunen-Foeve 
expansion.  [4,5] 

The  Karhunen-Foeve  expansion  provides  an  orthogonal  decomposition  for  any  second 
order  random  field  (stochastic  process)  w(x,a>)  where  x  is  the  space  coordinate  and  0  is  an 
element  of  the  random  event  space.  The  Karhunen-Foeve  expansion  for  this  random  process  is 
written  as 


oo 

W(x,6))  =  w(x)  +  ^n((0)j\fn(x)  (1) 

n= 0 

where  £,n ( co)  is  a  set  of  random  variables  defined  on  the  event  space  Q ,  and  f„(x)  is  a  set  of 
detenninistic  orthogonal  functions.  It  can  be  shown  that  the  orthogonal  function  set  consists  of 
the  eigenfunctions  for  the  covariance  function  of  the  random  process.  The  Xn  are  the  associated 
eigenvalues. [4,5]  The  space-dependent  mean  value  for  the  random  process  is  w(x).  As  defined, 
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this  expansion  is  convergent  in  mean  squared  error  and  is  unique,  but  it  is  mostly  of  theoretical 
interest.  It  is  not  very  useful  for  practical  computations  since  the  covariance  function  for  the 
random  process  must  be  known  a  priori  in  order  to  guarantee  the  expansion’s  convergence  and 
uniqueness. [4]  In  most  cases,  the  covariance  function  is  unknown.  Still,  the  expansion  motivates 
the  search  for  other  orthogonal  decompositions  that  are  more  easily  applied  even  if  the 
convergence  of  individual  decompositions  must  be  carefully  monitored.  The  method  showcased 
in  this  report  applies  a  polynomial  chaos  expansion  to  construct  the  orthogonal  decomposition. 
An  early  example  of  this  expansion  is  Weiner’s  homogeneous  chaos. [10] 

Wiener’s  continuous  homogeneous  chaos  cast  in  three  dimensions  is  a  measurable 
function  p  with 


P  =  P(xv  x2,x3;/?)  (2) 

In  this  chaos,  x-(xl,x2,xi)  is  the  detenninistic  space  point  while  /?  e  [0,1]  represents  the 

random  range  of  the  stochastic  process.  Other  chaos  fonnulations  have  been  proposed  and 
implemented.  A  specific  chaos  is  due  to  Cameron  and  Martin  and  is  denoted  as  a  Fourier- 
Hennite  chaos  for  a  chosen  functional.  The  Cameron-Martin  theorem  states  that  the  Fourier- 
Hennite  series  of  any  real  or  complex  functional  F[x]  in  L2  converges  to  F[x]  in  the  L2  sense 
with  Wiener  measure.  [3]  In  the  same  reference,  the  orthogonality  of  the  Fourier-Hennite  chaos  is 
also  demonstrated.  The  measure  (or  weighting  function)  used  in  integration  is  the  same  as  the 
probability  density  function  for  Gaussian  random  variables. [1 1]  This  chaos  is  still  widely  applied 
for  uncertainties  consisting  of  Gaussian  random  variables.  To  provide  for  the  assessment  of  other 
types  of  random  variables,  other  chaos  formulations  have  been  developed.  These  fonnulations 
are  included  in  the  generalized  polynomial  chaos  or  “Askey  Scheme”.[ll]  In  this  generalized 
chaos,  an  appropriate  set  of  orthogonal  polynomials  is  paired  with  a  random  variable 
distribution.  As  described  in  Reference  [1 1,12,13],  these  combinations  are  shown  in  Table  1. 

Table  1.  Distribution  Function/Orthogonal  Polynomial  Chaos  Combinations 


Distribution  Function 

Orthogonal  Polynomial  Set 

Gaussian 

Hermite 

Uniform 

Legendre 

Gamma 

Laguerre 

Beta 

Jacobi 

Poisson 

Charlier 

Negative  Binomial 

Meixner 

Binomial 

Krawtchouk 

Hypergeometric 

Hahn 

Consider  a  Hermite  chaos  for  random  function  a  formulated  for  exactly  one  random  variable  E, 
defined  for  the  event  space  Q  (where  co  e  Q  ).  For  this  scenario,  the  chaos  expansion  is  written 
as 
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(3) 


«(«»))  =  s  a,H,  ,«(<»)) 

n=  0 


The  form  of  this  expansion  is  quite  simple  since  there  are  no  space  or  time  parameters  explicitly 
associated  with  the  random  function.  In  this  sense,  this  problem  is  analogous  to  a  single  coin  toss 
or  dice  throw  problem.  Specifically,  the  an  are  expansion  coefficients  that  are  detennined  to 

complete  the  representation  of  the  random  function  a  .  The  notation  Hn  represents  the  order  n 
Hennite  polynomial  where  n  =  0, 1,  2  ... .  The  Hennite  polynomials  are  written  as  functions  of 
the  random  variable  £,  e.g.,  H0(%)  =  1;  HX{^)  =  E,\  H2(^)-^2- 1,  etc. [4]  If  the  random 
function  involves  more  than  one  random  variable,  e.g.,  E,x  and  £2,  then  the  expansion  functions 
are  formed  by  tensor  products  of  the  polynomials  cast  in  one  random  variable.  The  procedure  for 
determining  the  expansion  coefficients  an  is  explained  in  the  following  section  of  this  report. 
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3.0  METHODS,  ASSUMPTIONS  AND  PROCEDURES 


In  the  discussions  below,  we  develop  the  equations  needed  in  applying  a  polynomial 
chaos  to  represent  uncertainties  in  stochastic  processes.  The  two  principal  approaches  for 
applying  a  polynomial  chaos  are  (i)  the  intrusive  method  and  (ii)  the  non-intrusive  method. 
Although  this  report  concentrates  on  the  non-intrusive  application,  the  intrusive  approach  is  also 
briefly  described.  For  the  example  problems  shown  in  Section  3,  detailed  equations  are  presented 
for  the  Legendre  polynomial  chaos  associated  with  for  a  system  of  uniform  random  variables. 
Due  to  the  simplicity  of  the  probability  measure  for  the  uniform  random  variable,  this  chaos  is 
the  easiest  to  derive.  As  a  result,  the  development  presented  below  is  intended  to  be  pedagogical. 
As  a  second  discussion,  the  Hermite  polynomial  chaos  is  developed  for  Gaussian  random 
variables. 

3.1  Intrusive  Polynomial  Chaos 

This  manifestation  of  polynomial  chaos  is  denoted  intrusive  since  it  entails  alterations  to 
the  system  of  governing  equations.  An  example  discussed  in  the  literature  involves  solution  of 
the  incompressible  Navier-Stokes  equations. [1,1 7]  In  this  case,  a  single  parameter  is 
decomposed  say,  u,  the  x-velocity  component. 


u(x,t  =  0)  =  u(x)  +  %  u{x)  (4) 

In  this  expression,  %  is  a  random  variable  with  unit  variance.  Polynomial  chaos  expressions  for 
all  stochastic  quantities,  e.g., 


u(x,  t,  £)  =  JX(*, 

k= 0 


(5) 


Expressions  such  (5)  are  substituted  into  the  Navier-Stokes  equation.  By  use  of  the  Galerkin 
procedure,  evolution  equations  for  expansion  coefficients  uk  may  be  derived. [17]  For  instance, 


duk 

dt 


+ 


£2X-.(“,V>" 


i= 0  y'=0 


=  -Vpk  +  vV\ 


(6) 


for  k  =  0,...,P,  and 


V  •  uk  -  0 


Mijk  = 


ML 


(7) 


(8) 


The  Mjjk  are  convolution  integrals  rendering  numbers  for  use  in  (6).  Equations  (6)  and  (7)  must 
be  solved  via  numerical  means  to  obtain  the  expansion  coefficients.  This  alteration  of  the 
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governing  equations  is  characteristic  of  intrusive  polynomial  chaos  methods.  Given  the  difficulty 
associated  with  implementing  these  schemes,  implementation  details  are  not  provided  here. 

3.2  Fundamentals  of  Polynomial  Chaos 

As  its  most  basic  idea,  a  polynomial  chaos  expands  the  uncertainty  of  a  stochastic  process 
in  tenns  of  a  system  of  orthogonal  polynomials.  By  increasing  the  number  of  tenns  in  the 
expansion,  the  overall  accuracy  of  the  method  is  enhanced.  Of  course,  obtaining  the  expansion 
coefficients  for  the  series  constitutes  the  majority  of  the  labor  required  to  implement  this  method. 
The  polynomial  chaos  approach  applied  here  is  a  direct  application  of  the  equation 


a  =  (9) 

7=0 

where  a  (x,  /,  c )  is  the  random  process  of  interest;  VF/  is  an  element  of  an  orthogonal  family  of 

polynomials,  and  the  adx,  t)  are  the  expansion  coefficients  for  these  polynomials.  Of  course,  £, 

is  a  vector  of  random  variables  representing  the  system’s  uncertainties.  Each  of  these  random 
variables  is  usually  chosen  as  of  either  “standard”  unifonn  or  Gaussian  fonn.  Of  course,  other 
types  of  random  variables  can  be  used.  The  usual  procedure  is  to  have  the  random  variable  match 
the  random  characteristics  of  the  uncertainty.  vF/.(<|)is  a  multi-dimensional  orthogonal 
polynomial  in  that 


J'jy* (io) 

U(|) 

Equation  (10)  is  the  statement  of  orthogonality.  The  expected  value  for  a*  may  be  computed  as 
follows.  The  expected  value  is  defined  as  (axj ;  therefore,  by  truncating  the  series  to  a  finite 
number  of  terms, 


(a\x,  0)  =  YjOCj(x,  t )  J Tf(|)  w(|)  </| 

7=°  1 3(|) 


(ID 


(a\x,  0)  =  a0(x,  t )  J  TVf )  w(|)  d£,  +  $>,.(*,  t)  J  Tf(|)  w(|)  dl  (12) 

D(h  )  = 1  £>(f) 

Although  series  truncation  is  common  practice,  it  must  be  performed  carefully  since  it  may  affect 
series  convergence.  Equation  (12)  may  be  simplified  that  realizing  that  VF0  is  the  zeroth  order 
polynomial,  a  constant,  so  it  can  be  chosen  as  unity.  Thus, 
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a\x,  O)  =  a0(x,  t )  J  w(|)  d£  +  Z““AET“^  J  %  •  ^(1)  w(l)  <*l 

£>(|) 


VI/ 

/'=*  0  £>(f) 


(13) 


In  equation  (13),  the  first  tenn  is  evaluated  by  noting  that  the  probability  measure  is  written  as 


dw(£)  =  w(£) 


Applying  the  rules  of  probability, 


Jw(|)«?| 


=  1 


(14) 


(15) 


D(S) 


Next,  observe  that  the  second  term  contains  the  rewritten  integral 

|'V¥.(|)w(|)J|,  j  =  l,2,...,P 


(16) 


D(i) 


The  rewrite  is  made  possible  since  %  is  a  constant,  equal  unity.  By  using  (10)  and  (15),  the 
integral  in  (16)  is  zero.  Hence, 


a*(x,  t)^  =  a0(x,  t ) 


The  variance  can  be  computed  by  a  similar  set  of  arguments.  Recall  that 


var(or*)  =  (a*2) -(a* 


The  first  tenn  in  (18)  is  evaluated  by  substituting  (9),  i.e., 


a2)  =  i>/Fy(|)  •  5>***(l)  =  YL°j  a*^(l)^*(l) 


V=° 


k= 0 


\j= 0k=0 


By  using  the  integral  to  compute  the  expectation, 


(17) 


(18) 


(19) 


«’2) 

j=0  k=0  £,(1) 


Applying  (10),  we  obtain 


*2 


a 


p  p 


j=0  k-0 


7=0 


(20) 


(21) 
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Substituting  (21)  and  (17)  into  (18),  we  have  that 


var(«*  (x,  t))  =  a)  (x,  t)  ( Y? )  -  a]  (x,i t)  =  Y  a) (*>  0  (^j  ) 

j= 0  7=1 

The  standard  deviation  a(x,t)  is,  by  definition,  the  square  root  of  the  variance,  so 


(22) 


a(x,  t )  =  var(a*(x,  t ))  =  ^JY  aj  (*> 0  (xpy  }  (23) 

As  illustrated  by  (17)  and  (23),  the  mean  and  standard  deviation  for  the  random  process  are 
easily  calculated  from  the  polynomial  chaos  expansion. 

3.3  Point  Collocation  Non-Intrusive  Polynomial  Chaos 

Non-intrusive  polynomial  chaos  methods  do  not  require  modification  of  the  governing 
equations.  Of  equal  importance  is  the  fact  that  there  is  more  than  one  non-intrusive  method. 
Spectral  progression  and  linear  regression  are  two  such  methods.  [18]  As  with  the  intrusive 
method,  equation  (9)  is  applied  to  physical  properties  of  interest  say,  absolute  pressure  for  a  fluid 

dynamics  problem.  The  property  of  interest  is  represented  by  a*{x,t,^).  It  is  important  to 

realize  that  E,  represents  the  input  or  driving  uncertainty  say,  a  varying  initial  or  boundary 

condition.  For  the  point  collocation  method,  the  random  variable  E,  is  randomly  sampled  in 
accordance  with  its  statistical  distribution.  Each  sample  fonns  a  realization  of  the  random 
process  a* .  For  a  direct  detennination  of  the  polynomial  chaos  coefficients  a  ,  P+1  samples, 

,  j  -  0,1,..., P  ,  of  the  random  process  a*  are  computed  (perhaps  by  use  of  fluid  dynamics 
computer  codes  for  each  ^.).  Applying  this  idea  prior  to  series  truncation,  equation  (9)  is 
rewritten  as 


a\x,  t,  | ,.)  =  Yak(*>  0  ^(f,)  >  J  =  0,1,--- ,P 


(24) 


k= 0 


The  polynomials  yYk(Ej)  are  calculated  for  the  choice  of  Ef ,  so  these  parameters  are  known 
quantities.  Equation  (24)  actually  represents  a  system  of  P+1  equations  in  P+1  unknowns.  As  a 
system  of  linear  equations,  (24)  can  be  solved  for  the  ak .  The  example  problems  documented 
later  in  this  report  apply  LU-decomposition  with  partial  pivoting  to  solve  this  equation. [19] 
Naturally,  the  title  of  this  method  arises  from  the  use  of  the  collocation  points  E  • 

The  point  collocation  method  applied  here  is  introduced  in  References  14,  15,  20  and  21. 
It  is  relatively  easy  to  implement,  yet  it  is  subject  to  some  limitations  that  must  be  mentioned. 
Recall  from  the  preceding  paragraph  that  the  polynomial  chaos  equations  are  defined  at  the  P+1 
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collocation  points.  For  a  given  continuous  random  variable  distribution  (e.g.,  unifonn,  beta, 
Gaussian),  a  random  selection  of  collocation  points  is  not  unique.  For  this  reason,  the  computed 
mean,  standard  deviation  and  other  statistical  parameters  are  not  unique.  A  number  of  numerical 
techniques  have  been  designed  as  remedies  for  this  difficulty.  One  approach  is  to  generate  a 
more  even  point  distribution  to  provide  better  coverage  of  the  collocation  point  space.  This  may 
be  accomplished  through  the  use  of  Hammersley  or  Halton  data  sets. [22, 23]  Numerical 
experiments  conducted  by  the  author  demonstrate  that  a  Hammersley  point  set  does  provide 
better  coverage  of  the  uncertainty  distribution.  Unfortunately,  the  lack  of  true  randomness  in  this 
data  set  tends  to  induce  linear  dependence  within  the  system  (24).  The  system  determinant 
collapses  to  zero  causing  the  matrix  solution  to  fail.  Still,  this  difficulty  can  be  bypassed  by 
adding  a  pseudo-random  perturbation  to  the  Hammersley  points’  coordinates.  Computations 
have  shown  that  this  algorithm  works  reasonably  well. 

An  alternative  that  can  grant  good  coverage  of  uncertainty  space  (when  combined  with 
random  point  collocation)  is  the  method  of  least  squares.  In  this  approach,  a  superset  of  R+  1 
random  collocation  points  is  generated  where  R  >  P.  This  number  of  collocation  points  forces 
(24)  to  become  an  over  determined  system  of  equations.  The  over  determined  version  of  system 
(24)  may  be  written  as 

[¥]«  =  «*  (25) 

where  [lF]  is  an  (R  + 1)  x  (P  +  1)  matrix  of  the  fonn 

%(lo)  ^i(Io)  -  ^p(lo)  " 

%(li)  'J'i(li)  -  ¥,(!) 

m=  :  :  :  (26) 

%(!*-! )  -  VpiL i) 

%(4)  ^(4)  •••  %(4)_ 

and  a  is  a  (P  + 1)  x  1  vector  containing  the  unknown  expansion  coefficients, ,  i.e., 

a  =  (a0,  ax,  ...,ap)T  (27) 

The  right  hand  side  of  (25)  is  the  (R  + 1)  x  1  vector  denoted  a*  containing  the  values  of  the 
random  process  evaluated  at  the  R  +  1  collocation  points.  This  vector  is  written  as 

a  =(a*,  ax,...,aR)T  (28) 

The  classical  least  squares  system  of  equations  is  obtained  by  multiplying  (25)  by  the  transpose 
of  polynomial  matrix,  i.e., 

[Tf[T]a  =  [Tf  a  (29) 
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The  matrix  inner  product  [VP]7'[VP]  is  a  (P  +  1)  x  (P  +  1)  square  matrix.  Similarly,  [vF]r<5  is  a 
vector  with  P+ 1  elements.  With  these  dimensions,  (29)  is  solvable  by  standard  numerical  linear 
algebra  techniques.  The  specific  matrix-vector  forms  are 


mTm= 


yvj/2  yip  vp 

/  ,  1  Ok  Z^u  1  Ok  1  1 

k=0  k=0 

yxj/  xj/  yi}/2 

Z^  T\k  Tok  z^  Ti* 


k=0 

R 


k=0 


k=0 


Z'f'A  Z'1'.™'*'.. 


k=0 


[^fa*  = 


z*. 

k=0 

R 

Z'1'. 

k=0 

R 

Z'J', 

k= 0 


\kak 


Pkak 


iw 

k=0 

Z'f'u'P, 


k= 0 


k= 0 


2 

PA: 


(30) 


(31) 


The  classical  least  squares  formulation  works  well  when  applied  to  problems  involving  the 
uniform  probability  distribution.  One  may  recall  that  unifonn  probability  distribution  has  a 

constant  weight  function  w(d;)  within  its  probability  measure.  Other  probability  distributions 
have  non-constant  weight  functions.  An  examination  of  the  above  equations  reveals  that  the 
classical  least  squares  approach  does  not  incorporate  the  effect  of  the  weight  function  in  a  direct 
manner.  The  effect  of  a  non-uniform  distribution  is  conveyed  only  in  a  tacit  manner,  by  a  careful, 
histographical  selection  of  collocation  points.  To  provide  an  improved  capability  for  the  least 
squares  approach,  an  alternative  fonnulation  does  include  the  effects  of  weighting  the  collocation 
points. 


Rigorous  least-squares  derivations  are  based  upon  a  minimization  of  squared  error.  For 
the  point  collocation  polynomial  chaos  problem,  the  error  may  be  formulated  as  follows. [24,25] 
Recall  that  the  polynomial  chaos  representation  for  random  process  a  is  written  as 

i=0 

For  this  analysis,  the  truncated  series  is  employed.  The  total  squared  error  E  for  the  chaos 
expansion  may  be  expressed  as  a  function  of  the  expansion  coefficients  a  ,  j  =  0, 1, . . .,  P ,  i.e., 
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E{a0,a„...,aP)=  J  a*(|) 


Equation  (33)  integrates  over  i;  as  a  continuum.  To  evaluate  the  squared  error  for  the  collocation 
points,  we  apply  the  Monte  Carlo  average. [24,25] 


E(a„a„...,P )  =  -t-Xf  °'(l>  -  )  »<l, ) 

TV  +  1 


To  minimize  the  squared  error,  we  set  a  condition  on  the  first  partial  derivative. 


=  0,  k  =  0,\,...,P 


By  applying  (35)  to  (34),  we  obtain 


r)F  1  Q  f  p  \ 

—  =  7 a(l)-£a,VXl)  w(l)  =  0,  *  =  0,1,...„P 
dak  Q  dak  \  to  j 


After  differentiation  and  substitution,  we  exchange  the  order  of  summation;  the  result  is 


;=o  W= o 


2>,  Xiw#, )▼,(#, Ml,)  -£«'<#, wl, Ml, ),  * = o.i,. ...p 


Equation  (37)  constitutes  a  system  of  linear  equations.  By  setting 


the  linear  system  may  be  written  as 


*-=o,i . p 


Clearly,  the  probability  weight  function  is  incorporated  in  both  (38)  and  (39).  The  system  (40) 
can  be  solved,  like  (29),  by  methods  such  as  LU-decomposition. 
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3.4  Legendre  Polynomial  Chaos 


For  this  chaos,  the  uncertainties  are  represented  by  a  random  vector  E,  comprised  of  a  set 
of  n  random  variables,  i.e., 


|  =  (£,&,.  (41) 

In  addition  to  the  random  character  contained  in  (4),  the  random  process  a  is  also  defined  in 
space  and  time. [14,  15]  This  process  is  written  as 

a(x,  t,  g)  =  '£iaj(x,t)'¥j(£)  (42) 

j= o 

Note  that  this  polynomial  decomposition,  in  full  fonn,  possesses  an  infinite  number  of  terms  like 
a  typical  Fourier  series.  For  practical  computations,  the  series  is  truncated  as  in  (42)  to  retain 
only  a  finite  number  of  terms.  The  maximum  number  of  terms,  P  +  1 ,  may  be  computed  from  the 
formula 


P  +  1 


(P  +  n)l 

p\  n\ 


(43) 


where  p  is  the  order  of  the  Legendre  polynomial  set  used  for  the  chaos.  This  total  number  of 
expansion  terms  can  be  arrived  at  by  combinatorial  arguments.  Therefore,  if  p  equals  five,  then 
Legendre  polynomials  of  orders  zero  through  five  are  used  to  form  the  chaos.  The  construction 
of  expansion  (42)  is  very  important;  each  tenn  is  separated  into  a  detenninistic  coefficient 
0Cj(x,  t)  and  a  random  tensor  product  polynomial  ++(2 ) .  The  tensor  product  polynomial  is  a 
product  of  single  random  variable  Legendre  polynomials.  That  is  to  say, 


(I )  =  Lpl(& )Lpl  (&)•••  Lpn  (^  ) 


(44) 


Table  2.  Legendre  Polynomials 


Order 

Legendre  Polynomial 

0 

A>(£)  = 1 

1 

A(£)  =  £ 

2 

4(£)  =  -i 

3 

4(£)  =  £(5£2-3) 

4 

La{%)  =  3-30<f +35£4 

5 

4(£)  =  £(15-70£2  +  63£4) 
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Legendre  polynomial  orders  pi,  p2,...,  pk  are  constrained  such  that  0  <  pk  <  p ,  1  <  k  <  n  .  An 
individual  Legendre  polynomial  Lpk  has  algebraic  order  k  .  The  first  six  Legendre  polynomials 

are  presented  in  Table  2.  Of  course,  these  polynomials  are  unique  up  to  a  constant  multiple.  With 
some  effort,  higher  order  Legendre  polynomials  may  be  calculated  by  using  the  Rodrigues 
fonnula 


4(£>  = 


1  rT(cf-l)" 
2"  n\  d£n 


(45) 


It  follows  from  (44)  that  TV  has  the  algebraic  order  p\  +  p2  +  ■■■  +  pk.  It  is  computationally 
advantageous  to  arrange  the  terms  in  (42)  by  increasing  algebraic  order  beginning  with  the  zeroth 
order  term,  then  the  first  order  terms  and  so  forth.  For  two  random  variables  4\  and  g2 ,  this 
ordering  of  the  VF/  can  be  determined  from  organization  shown  in  Table  3.  Those  VF/  possessing 


Table  3.  Organization  of  Legendre  Polynomials  for  Two  Random  Variables 


%=l0(4)l0(4) 

%=L0(4)L1(4) 

w5=l0(4)l2(4) 

V  =  L0(4)LP(4) 

'4  =  44)44) 

^=4(^)4^) 

T  =  4(^(^) 

4  =4(^)4^) 

%=  4(^)4  (4) 

'F  =  4(4>4(4> 

'¥  =  l2(4)lp(4) 

'¥  =  Lp(4)L0(4) 

'F=z,(£)4(£) 

'¥  =  lp(4)l2(4) 

=lp{4)lp{4) 

the  same  algebraic  order  are  located  on  diagonals  (lower  left  to  upper  right)  in  this  table.  The 
first  diagonal  has  order  zero;  the  second  diagonal  ('Fj  and  4f  )  has  order  one;  the  third  diagonal  ( 
Tf ,  T4  and  Tf )  has  order  two,  and  the  pth  diagonal  has  order  p- 1.  The  remaining  diagonals 
follow  suit,  culminating  with  the  bottom  right  comer  table  entry  with  order  p2  - 1 .  The 
numbering  of  the  'P  follows  that  shown  in  the  Table  3,  but  the  larger  j  indices  are  omitted 

because  of  the  lengths  of  the  expressions.  Still,  individual  indices  are  not  difficult  to  discern 
because  the  length  of  each  diagonal  is  known.  That  is,  the  number  of  tenns  possessing  the  same 
algebraic  order  is  the  same  as  the  number  of  elements  in  the  associated  diagonal. 

To  support  accurate  calculations  of  statistical  properties  via  equations  (17),  (22)  and  (23), 
it  is  necessary  to  normalize  the  Legendre  polynomials.  The  nonnalized  magnitude  “norm”  for  an 
orthogonal  polynomial  / (£)  may  be  calculated  as  follows.  The  squared  norm  of  / (£)  ,  denoted 

(y2) ,  is  given  by 


(f2)=  J/2(<T)<Mf) 

O(f) 


(46) 
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In  equation  (46),  Q(<^)  is  the  domain  of  the  random  variable  £  while  dw(%)  is  the  probability 
measure  associated  with  the  random  variable.  The  Legendre  chaos  is  associated  with  a  standard 
uniform  random  variable  2  where  d,  e  [-1, 1] .  The  mean  of  this  random  variable  is  zero,  and  the 
probability  measure  is  given  by 


dw(Z) 


(47) 


With  the  use  of  (47),  the  squared  nonns  for  the  Legendre  polynomials  found  in  Table  2  are 
computed  and  recorded  in  Table  4. 


Table  4.  Squared  Norms  for  the  Legendre  Polynomials  in  One  Random  Variable 


Order 

Squared 

Norm 

0 

1 

1 

1/3 

2 

4/5 

3 

4/7 

4 

64/9 

5 

64/11 

3.5  Hermite  Polynomial  Chaos 

A  distinct  advantage  of  the  polynomial  chaos  method  is  the  generality  of  its  mathematical 
form.  The  equation  for  the  expansion,  the  structure  of  the  uncertainty  vector  and  number  of  terms 
in  the  expansion  remain  unchanged.  Only  the  orthogonal  polynomial  set  and  the  probability 
weight  function  (measure)  change.  The  Hennite  chaos  is  based  upon  the  standard  nonnal 
(Gaussian)  distribution,  i.e., 


"<l)  =  7rLTexp(i|.|)  (48) 

(In) 

Equation  (48)  is  suited  for  an  uncertainty  vector  of  arbitrary  integer  length  n.  For  a  chaos 
involving  only  one  random  variable  £  (the  uncertainty), 


"<£)  = 


exp 


(  £2  \ 


(49) 


Similarly,  for  two  random  variables  £  =  (q, ,  42 ) ,  the  weight  function  is 
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(50) 


Hi) 


f  zl 


2  n 


exp 


2\ 


The  Hennite  polynomials  constitute  the  orthogonal  spanning  set  required  for  this  chaos.  The 
Hermite  polynomials  for  one  random  variable  may  be  calculated  from  the  Rodrigues  fonnula 


(1 

■  dj 

(  1  Al 

-1  y— 

exp 

\2  j 

d£J 

K  2  )_ 

(51) 


The  first  five  of  these  polynomials  are  listed  in  Table  5  along  with  their  squared  norms. 


Table  5.  Hermite  Polynomials 


Order 

Hermite  Polynomial 

Norm 

0 

nm  =  1 

1 

1 

HM)  =  S 

1 

2 

H^)=e- 1 

2 

3 

6 

4 

H<(£)  =  3-6<f +r 

24 

For  systems  involving  more  than  one  Gaussian  random  variable,  this  chaos  is  still  applicable,  but 
the  calculation  of  Hennite  polynomials  is  more  difficult.  These  computations  are  reserved  for  a 
future  report.  The  chaos  formulations  presented  in  this  section  and  the  previous  one  are  well 
suited  for  the  example  problems  that  follow. 
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4.0  RESULTS  AND  DISCUSSION 


In  this  section,  four  test  problems  are  described  for  the  polynomial  chaos  solution.  The 
first  three  problems  are  designed  for  solution  by  the  Legendre  chaos.  The  final  problem  applies 
the  Hermite  chaos  to  resolve  the  uncertainty  associated  with  shocked  gas  dynamics 
configuration.  The  exact  statistical  parameters  for  the  first  four  cases  are  derived  in  the 
appendices. 

4.1  Test  Problem  1 

This  test  case,  designed  by  the  author,  involves  a  random  process  operating  on  one 
random  variable  %  e  [-1, 1] .  The  random  process  / (£)  is 


/(£>  =  i-i(r  +  r  +  f) 


(52) 


Statistical  parameters  for  this  random  process  are  computed  through  the  use  of  a  fifth  order 
Legendre  polynomial  chaos.  The  exact  solution  for  this  problem  is  provided  in  Appendix  A.  In 
the  terminology  of  Section  2,  p  =  5,  and  n  =  1 ,  so 


P  + 1 


5!  1! 


(53) 


For  a  direct  polynomial  chaos  solution,  a  system  of  6  equations  in  6  unknowns  must  be  solved. 
The  equations  are  formulated  at  the  6  collocation  points  shown  in  Table  6.  The  mean  and 
standard  deviation  calculated  for  this  test  problem  are  shown  in  Table  7  and  compared  with  the 
exact  calculations.  The  agreement  between  the  chaos  estimates  and  the  exact  parameters  is  quite 
good  especially  given  the  low  number  of  collocation  points. 

Table  6.  Collocation  Point  Coordinates  for  Test  Problem  1 


Point  No. 

4 

1 

-0.90 

2 

-0.60 

3 

-0.30 

4 

0.20 

5 

0.55 

6 

0.95 

Table  7.  Statistical  Parameters  for  Test  Problem  1 


Parameter 

5th  Order  Chaos 

Exact 

Mean 

0.77238 

0.7746 

Standard  Deviation 

0.26228 

0.26258 
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4.2  Test  Problem  2 


This  test  problem,  again  developed  by  the  author,  is  cast  in  terms  of  a  random  process  / 
operating  on  one  random  variable  e  [-1,1]  ■  Specifically, 


/«)  =  t(l-f2)(sin(10f)  +  3) 


(54) 


The  fifth  order  Legendre  chaos  is  also  applied  in  estimating  the  mean  and  standard  deviation  for 
this  random  process.  As  shown  for  the  previous  case,  six  collocation  points  are  required  for  this 
test  problem.  The  same  list  of  collocation  points  provided  in  Table  6  is  used  for  this  problem. 
The  Legendre  chaos  is  solved  for  this  test  case,  and  the  estimated  statistical  parameters  are  listed 
and  compared  against  the  exact  parameters  in  Table  8. 

Table  8.  Collocation  Point  Coordinates  for  Test  Problem  2 


Parameter 

5th  Order  Chaos 

Exact 

Mean 

0.49401 

0.50 

Standard  Deviation 

0.34134 

0.26048 

For  the  low  number  of  collocation  points  applied  here,  the  agreement  between  the  estimated  and 
exact  parameters  is  reasonably  good. 

4.3  Test  Problem  3 

The  first  two  test  problems  addressed  by  this  report  involve  one  random  variable  and 
serve  to  verify  the  polynomial  chaos  solution  algorithms  discussed  in  Section  2.  This  test 
problem  is  cast  in  two  random  variables  and  has  a  mixture  of  algebraic  and  transcendental 
functions. [15]  Specifically, 


f(xl,x2)  =  ln(  1  +  x2 )  sin(5x2 )  (55) 

where  jq  and  x1  are  uniform  random  variables  with  mean  2.0  and  a  coefficient  of  variation 
(CoV)  of  20%.  The  mean  and  CoV  can  be  used  to  identify  the  domain  of  x,  and  x2 ,  i.e., 

Xj,x2  e  [a,b]  =  [2.0 - 0.4V3,2.0  +  0.4-V3]  (56) 


Random  variables  xl  and  x2  are  regarded  as  functions  of  standard  uniform  random  variables  ^  2 

with  mean  zero  and  standard  deviation  V3  .  That  is,  £e[-l,l].  The  domain  of  <^12  can  be 
mapped  onto  the  domain  of  Xj  2  (and  vice  versa)  by  the  transformation 


x  ■=  a  + 


(. b-a ) 


(57) 
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The  exact  statistical  parameters  for  this  random  process  are  described  in  Appendix  C.  The 
number  of  terms  required  for  the  Legendre  polynomial  chaos  is  calculated  as 

p  +  l  =  (5  +  2)!  =2i  (58) 

5!2! 

This  polynomial  chaos  is  solved  for  two  different  point  distribution  schemes.  The  first  method  of 
distributing  collocation  points  is  the  random  method.  A  linear  congruential  random  number 
generator  is  used  to  distribute  points  in  the  interval  [0,1]. [26]  Then  these  points  are  mapped  onto 
the  domain  of  %  for  use  in  the  random  process  formulas  (55)  and  (57).  A  direct  solution  is 
perfonned  for  21  collocation  points.  These  points  are  listed  in  Table  9.  Secondly,  a  least  squares 
solution  is  accomplished  for  100  collocation  points.  For  brevity,  these  points  are  not  listed  here. 
The  fifth  order  direct  and  least  squares  solutions  for  this  problem  are  listed  and  compared  with 
exact  values  in  Table  10. 

Table  9.  Random  Collocation  Points  for  the  Direct  Solution  of  Problem  3 


Point 

£ 

4 

Point 

6 

6 

1 

0.9999974 

0.9555932 

12 

0.4861153 

0.1398800 

2 

0.6545904 

-0.2987570 

13 

0.9636399 

-0.1042775 

3 

0.7917758 

-0.6235767 

14 

-0.5915973 

-0.9762128 

4 

-0.4539059 

-0.7957101 

15 

0.7914310 

-0.4188466 

5 

0.5000229 

-0.1146452 

16 

0.4447752 

-0.6623921 

6 

-0.8423904 

-0.0558344 

17 

-0.8242336 

-0.8933002 

7 

-0.4100209 

0.7791101 

18 

0.3031493 

-0.9700800 

8 

0.5027170 

-0.8346537 

19 

-0.1349102 

0.5640060 

9 

-0.0248729 

-0.0391431 

20 

-0.7509418 

0.9213717 

10 

0.1213641 

-0.2343712 

21 

-0.5064132 

0.7136617 

11 

0.9234417 

0.2839582 

Table  10.  Statistical  Parameters  for  Test  Problem  3  with  Random  Point  Distributions 


Parameter 

Direct 

Least  Squares 

Exact 

Mean 

0.0729048 

0.078819 

0.079169 

Standard  Deviation 

1.14888 

1.11923 

1.12413 

A  second  method  of  solving  this  problem  involves  the  use  of  Hammersley  sampling  points.  Both 
the  direct  and  least  squares  algorithms  are  applied  in  estimate  the  statistical  parameters  for  this 
case.  For  a  direct  linear  solution,  the  21  Hammersley  collocation  points  are  listed  in  Table  11. 
Although  not  listed  here,  100  Hammersley  points  are  used  for  a  least  squares  solution.  Estimates 
of  the  mean  and  standard  deviation  for  problem  3  based  upon  the  Hammersley  point  distributions 
are  compared  against  the  exact  values  in  Table  12. 
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It  is  interesting  to  note  that,  on  the  whole,  polynomial  chaos  estimates  are  reasonably 
accurate  for  this  problem,  but  the  accuracy  of  the  least  squares  estimates  for  random  point 
distributions  compare  more  favorably  to  the  exact  values.  It  is  a  little  shocking  to  see  that 
estimates  based  upon  the  Hammersley  distributions  perfonn  more  poorly,  even  though  the 
Hammersley  distribution  provides  better  coverage  of  the  sample  space.  The  least  squares 
analysis  improves  the  estimates  made  by  using  both  random  and  Hammersley  point  distributions. 

Table  11.  Hammersley  Collocation  Points  for  the  Direct  Solution  of  Problem  3 


Point 

£ 

6 

Point 

£ 

6 

1 

0.9999974 

0.9555932 

12 

0.4861153 

0.1398800 

2 

0.6545904 

-0.2987570 

13 

0.9636399 

-0.1042775 

3 

0.7917758 

-0.6235767 

14 

-0.5915973 

-0.9762128 

4 

-0.4539059 

-0.7957101 

15 

0.7914310 

-0.4188466 

5 

0.5000229 

-0.1146452 

16 

0.4447752 

-0.6623921 

6 

-0.8423904 

-0.5583447 

17 

-0.8242336 

-0.8933002 

7 

-0.4100209 

0.7791101 

18 

0.3031493 

-0.9700800 

8 

0.5027170 

-0.8346537 

19 

-0.9700800 

0.5640060 

9 

-0.0248729 

-0.0391431 

20 

-0.7509418 

0.9213717 

10 

0.1213641 

-0.2343712 

21 

-0.5064132 

0.7136617 

11 

0.9234417 

0.2839582 

Table  12.  Statistical  Parameters  for  Test  Problem  3  with  Hammersley  Point  Distributions 


Parameter 

Direct 

Least  Squares 

Exact 

Mean 

0.061475 

0.069076 

0.079169 

Standard  Deviation 

1.183351 

1.123245 

1.12413 

4.4  Test  Problem  4 

This  test  case  is  the  last  of  the  truly  “academic”  demonstration  problems  addressed  here. 
The  random  process  to  be  considered  operates  on  four  independent  uniform  random  variables. 
The  random  process  is  defined  as 

f(x1,x1,xi,x4)  =  exp[l  ,5(xj  +x2  +  x3  +  x4 )]  (59) 

where  random  variables  xt  =  xi(^i) ,  i  =  1,  ...,  4  are  functions  of  the  standard  uniform  random 
variables  £ ,  where  f  e  [-1,1]  .[15]  Random  variables  xi  have  a  mean  value  of  0.4  and  a 
coefficient  of  variation  (CoV)  of  0.4.  With  the  use  of  this  infonnation,  it  can  be  shown  that  the 
domain  for  xi  is 


[a,b]  =  [0.4  -  0.08^12,0.4  +  0.08^12] 


(60) 


20 

Distribution  A 


For  a  fifth  order  polynomial  chaos  solution,  the  number  of  terms  for  the  expansion  is  calculated 
as  follows. 

(5  +  4)' 

P  +  1  =  ± - —  =  126  (4-10) 

5!4! 

A  tenn  in  the  expansion  polynomial,  following  equation  (44),  has  the  form 


=  4i(6)  LAQ  K 3(6)  44(6)  (4- 10a) 

where  0  <  nj  <  5 ,  j  =  1, . . . ,  4 .  In  this  case,  a  large  number  of  tenns  is  required  to  determine  the 
statistical  parameters  with  fifth  order  accuracy.  This  outcome  is  referred  to  as  the  “curse  of 
dimensionality”. [2 7]  That  is,  for  a  given  level  of  accuracy,  the  number  of  expansion  terms 
increases  significantly  for  each  additional  random  variable.  As  with  the  previous  test  cases,  the 
fifth  order  Legendre  chaos  is  applied  to  estimate  the  mean  and  standard  deviation  for  equation 
(59).  The  collocation  points  are  too  numerous  to  list  here.  Both  the  direct  solution  (126 
collocation  points)  and  the  least  squares  solution  (200  collocation  points)  have  been  computed 
for  this  problem.  The  results  are  presented  in  Table  13.  Because  of  the  poorer  performance  of  the 
Hammersley  point  set  used  in  the  previous  test  case,  it  is  not  used  here.  For  the  random  point 
sets,  the  results  are  quite  good  with  the  least  squares  case  granting  only  minor  improvements 
over  the  direct  method.  The  results  for  these  four  test  problems  indicate  that  the  Legendre  chaos 
performs  well  by  rendering  sound  estimates  for  statistical  parameters. 

Table  13.  Statistical  Parameters  for  Test  Problem  4  with  Random  Point  Distributions 


Parameter 

Direct 

Least  Squares 

Exact 

Mean 

12.3519 

12.3602 

12.3609 

Standard  Deviation 

6.1631 

6.1569 

6.1556 

4.5  Test  Problem  5 

The  preceding  problems  test  the  basic  operation  of  the  polynomial  chaos  algorithms 
given  known,  detenninistic  function  forms.  The  uncertainty  in  the  system  is  provided  by  the 
random  variable(s)  providing  input  data.  This  test  case  is  an  actual  application  of  this  idea.  Here, 
the  input  uncertainty  is  provided  by  changing  a  physical  boundary  condition  for  a  computational 
fluid  dynamics  (CFD)  solution. [14,  15,  21]  The  algebraic  functions  used  in  the  previous  test 
cases  are  replaced  by  separate  CFD  solutions  for  each  “randomly”  chosen  boundary  condition. 
The  Hermite  chaos  is  applied  to  solve  the  classical  oblique  shock  wave  problem.  A  Mach  3  flow 
field  of  a  calorically  perfect  gas  encounters  a  wedge.  To  satisfy  this  solid  surface  inviscid 
boundary  condition,  the  flow  must  turn  at  the  cusp  of  the  wedge  or  ramp.  Physics  dictates  that 
the  turning  occurs  through  an  oblique  shock  wave  originating  the  wedge’s  cusp.  This  scenario  is 
illustrated  in  Figure  1.  The  uncertainty  is  provided  by  perturbations  in  the  wedge  angle.  The 
mean  wedge  angle  0  is  five  degrees,  and  the  distribution  of  possible  wedge  angles  is  Gaussian 
with  a  coefficient  of  variation  of  10%. [21]  Based  upon  this  infonnation,  the  standard  deviation 
for  the  wedge  angle  distribution  is 
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M-  Mach  3 


<j  =  0.100  =  0.5° 

The  probability  distribution  for  the  wedge  angle  is  given  by 


(63) 


w(0)  = — ^=exp 
(jyJ2  n 


(0-0? 

2<j2 


Hence, 


1  00 


rV2 ~n 


00 

f  „  „ 

{o-of  1 

J  exp 

—  00 

1 

(N 

b 

(N 

1 

_ i 

dO  = 1 


(64) 


(65) 


By  using  the  transformation 


£  = 


e-e 


(66) 


we  obtain  the  standard  nonnal  distribution 

1  f 

^rJexp 


^2k  _ 


f  ^2\ 

v’^y 


d%  =  1  =>  w(£)  =  -7==exp 

^\2n 


f  gi\ 

v 


Equation  (66)  can  be  rewritten  as  follows. 


(67) 


0  =  0 +%c 7 


(68) 
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By  comparing  (64)  and  (67),  we  note  that  the  standard  deviation  for  the  standard  normal 
distribution  is  unity,  so  to  properly  sample  wedge  angles  for  this  problem,  samples  f  ■  must  be 

chosen  so  that  -2  <  <  2.  A  fourth  order  Hennite  chaos  may  be  used  to  estimate  statistical 

properties  for  pressure  throughout  the  flow  field.  Since  only  one  random  input  variable  (wedge 
angle)  is  involved,  the  total  number  of  expansion  terms  is  computed  as  follows. [21] 


P  + 1  = 


(4  +  1)! 
4!  1! 


=  5 


(69) 


With  the  use  of  these  equations,  the  sampled  wedge  angles  are  listed  in  Table  13. 

Table  14.  Sampled  Wedge  Angles  for  the  Oblique  Shock  Wave  Case 


i 

Sj 

Oj 

i 

-2 

4.0 

2 

-1 

4.5 

3 

0 

5.0 

4 

1 

5.5 

5 

2 

6.0 

To  statistically  sense  how  the  flow  field  is  perturbed  by  the  boundary  fluctuation,  the  flow  field 
is  sampled  at  the  three  points  listed  in  Table  14.  Points  1  and  3  differ  for  each  wedge  angle  since 
point  1  is  situated  in  the  shock  wave  while  point  3  is  on  the  wedge  surface.  As  a  result,  the  y 
coordinate  of  this  point  differs  for  each  wedge  angle.  In  Reference  15,  a  fixed  set  of  sampling 
points,  a  total  of  three,  is  used  for  all  wedge  angles.  This  approach  is  not  applied  here  since  an 
undisclosed  interpolation  algorithm  is  applied  in  Reference  15  to  refocus  the  data  at  the  sample 
point  locations.  Since  the  points  in  Table  14  are  relatively  close  to  the  points  used  in  the  archival 
reference,  we  assume  that  their  values  are  representative  of  the  random  process  at  the  fixed 
points. 


Table  15.  Sample  Point  Locations 


Wedge  Angle  (°) 

4.0 

4.5 

5.0 

5.5 

6.0 

Point 

X 

y 

y 

y 

y 

y 

1 

0.8984 

0.3669 

0.3675 

0.3822 

0.3879 

0.3969 

2 

0.8984 

0.2321 

0.2333 

0.2312 

0.2327 

0.2311 

3 

0.8984 

0.0644 

0.0722 

0.0813 

0.0879 

0.0985 

The  CFD  solution  for  each  wedge  angle  is  computed  by  using  the  Large  Eddy  Simulation 
with  Linear  Eddy  modeling  in  3  Dimensions  (LESLIE3D)  multiphase  physics  computer 
program.  LESLIE3D  is  developed  by  Professor  Suresh  Menon  at  the  Georgia  Institute  of 
Technology.  It  is  a  multi-block,  massively  parallel  computer  program  with  extensive  physics 
capabilities.  Structured  grids  are  developed  for  each  wedge  angle,  and  numerical  solutions  are 
computed  for  freestream  Mach  number  3.0,  sea  level  atmospheric  pressure  and  temperature 
300°K.  A  plot  of  the  velocity  field  for  the  5°  wedge  is  shown  in  Figure  2.  LESLIE3D  resolves  the 
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Figure  2.  Oblique  Shock  Solution  X-Velocity  Plot  for  the  5°  Wedge 


shock  wave  accurately  confirming  the  shock  angle  mandated  by  the  0  —  J3  —  M  relationship. [27] 
For  the  5°  wedge,  the  Cartesian  x-velocity  field  is  show  in  Figure  2.  The  pressure  field  is 
recorded  at  the  sample  points  (listed  in  Table  14  and  shown  in  Figure  2)  then  post-processed  by 
using  the  Hennite  polynomial  chaos.  The  ratio  of  pressure  versus  freestream  pressure  is  the 
random  process  considered  by  the  chaos.  The  computed  mean  and  standard  deviation  are  shown 
in  Tables  15  and  16,  respectively,  for  this  problem  and  compared  with  archived  results  taken 
from  Reference  15. 


Table  16.  Means  Computed  for  PI Pref  via  4th  Order  Hermite  Chaos 


Sample  Point 

Mean 

Archived 

1 

1.2091 

1.11433 

2 

1.4868 

1.45484 

3 

1.4651 

1.45472 

Table  17.  Standard  Deviations  Computed  for  P /  Pref  via  4th  Order  Hermite  Chaos 


Sample  Point 

Standard 

Deviation 

Archived 

1 

0.07435 

0.14033 

2 

0.07036 

0.05242 

3 

0.09073 

0.05227 

Overall,  the  agreement  between  the  present  Hermite  chaos  estimates  and  the  archived  estimates 
is  good  but  compares  less  favorably  than  the  other  test  problems.  There  are  a  number  of  reasons 
for  this  less  than  inspirational  comparison.  Our  method  utilizes  a  different  CFD  computer 
program  and  grid  than  does  the  archival  reference.  Also,  the  LESLIE3D  solution  is  perfonned 


24 

Distribution  A 


for  an  Oxygen-Nitrogen  mixture  with  air  mass  fractions  and  a  ratio  of  specific  heats  of  1.35. 
Also,  the  archival  solution  employs  interpolation  to  fixed  sample  points.  This  interpolation 
scheme  is  not  documented,  so  it  cannot  be  replicated  for  the  present  calculations.  For  future 
problems,  the  author  will  devise,  implement  and  document  an  interpolation  procedure  to  bridge 
the  gap  between  different  stochastic  realizations  of  a  more  complicated  problem  of  choice. 
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5.0  CONCLUSIONS 


This  expository  report  has  described  basic  aspects  of  the  polynomial  chaos  method  for 
representing  stochastic  or  random  processes.  These  methods  evolve  from  the  homogeneous 
chaos  introduced  by  Wiener.[10]  In  summary,  random  processes  can  be  represented  by  finite 
series  of  orthogonal  polynomials.  The  bulk  of  this  work  concentrates  on  the  point  collocation 
non-intrusive  polynomial  chaos  method  due  to  researchers  such  as  Walters  and  Hosder.[15]  Non- 
intrusive  methods  permit  statistical  estimates  without  altering  the  governing  equations.  Instead, 
polynomial  expansions  are  developed  at  a  set  of  collocation  points  selected  for  a  set  of  input 
random  variables.  These  random  variables,  representing  the  system  uncertainties,  drive  the 
random  aspects  of  the  overall  random  or  stochastic  process.  Although  the  point  collocation 
method  renders  non-unique  results,  its  estimates  have  shown  to  be  sufficiently  accurate  for  the 
included  test  problems.  In  each  case,  the  present  calculations  agree  well  with  archived  (or  exact) 
estimates. 

Polynomial  chaos  methods  are  not  without  difficulties.  Perhaps  the  most  prominent 
difficulty  is  that  of  dimensionality.  The  combination  of  polynomial  order  and  the  number  of 
random  variables  (uncertainties)  requires  an  increasing  number  of  expansion  terms.  To  resolve 
the  coefficients  for  a  system  of  n  terms,  an  n  x  n  matrix  must  be  solved.  Matrix  solutions  do 
require  significant  computer  resources  and  time.  Moreover,  for  a  time  and  space  dependent 
system,  the  chaos  must  be  solved  at  each  space  location  and  time  of  interest.  It  is  in  this  sense 
that  the  polynomial  chaos  method  begins  to  approach  the  workload  involved  in  a  Monte  Carlo 
simulation.  In  addition,  the  point  collocation,  non-intrusive  polynomial  chaos  method  produces 
results  that  are  not  unique.  Rather,  the  results  depend,  somewhat,  on  the  choice  of  collocation 
points.  The  analyst  must  ensure  that  the  random  space  receives  sufficient  coverage.  Moreover,  it 
must  possess  the  appropriate  level  of  randomness  with  the  appropriate  distribution  for  the 
uncertainties.  Other  analyses  have  shown  that  these  methods  have  difficulties  in  addressing 
higher  magnitude  uncertainties,  and  there  are  questions  concerning  the  accuracy  of  polynomial 
expansions. [28]  In  other  instances,  the  convergence  of  polynomial  chaos  expansions  may  be 
slow,  and  a  given  expansion’s  accuracy  may  not  increase  with  the  inclusion  of  more  terms. 
Statistical  moments  of  third  or  higher  order  may  not  be  accurate.  [29] 

In  a  future  report,  these  investigations  will  continue.  The  non-intrusive  polynomial  chaos 
method  will  be  applied  to  a  series  of  gas  dynamics  problems  involving  strongly  shocked  flow 
fields.  The  problems,  their  geometries  and  associated  physics  will  have  greater  complexity  that 
any  of  the  test  cases  shown  in  the  present  work.  A  number  of  uncertainties  will  be  examined, 
particularly  those  involving  the  boundary. 


26 

Distribution  A 


REFERENCES 


1.  Najm,  H.N.,  “Uncertainty  quantification  and  polynomial  chaos  techniques  in  computational 
fluid  dynamics”,  Annual  Review  of  Fluid  Mechanics,  Vol.  41,  2009,  pp.  35-52. 

2.  Ross,  S.M.,  Introduction  to  Probability  Models,  4th  Ed.,  Academic  Press,  Boston, 
Massachusetts,  1989. 

3.  Cameron,  R.H.  and  Martin,  W.T.,  “The  orthogonal  development  of  non-linear  functionals  in 
series  of  Fourier-Hermite  functionals”,  Annals  of  Mathematics,  Vol.  48,  No.  2,  pp.  385-392. 

4.  Ghanem,  R.G.  and  Spanos,  P.D.,  Stochastic  Finite  Elements,  A  Spectral  Approach,  Rev.  Ed., 
Dover  Publications,  Mineola,  New  York,  2003. 

5.  Pettersson,  M.P.,  Uncertainty  Quantification  and  Numerical  Methods  for  Conservation  Laws, 
Doctoral  Dissertation,  Institute  for  Computational  and  Mathematical  Engineering,  Stanford 
University,  2012. 

6.  Hammersley,  J.M.  and  Handscomb,  D.C.,  Monte  Carlo  Methods,  Methuen’s  Monographs  on 
Applied  Probability  and  Statistics,  Bartlett,  M.S.,  Ed.,  Fletcher  &  Son,  Ltd.,  Norwich,  Great 
Britain,  1964. 

7.  Anderson,  J.B.  and  Long,  L.N.,  “Direct  Monte  Carlo  simulation  of  chemical  reaction  systems: 
prediction  of  ultrafast  detonations”,  Journal  of  Chemical  Physics,  Vol.  118,  No.  7,  2003,  pp. 
3102-3110. 

8.  Johansen,  A.M.  and  Evers,  L.,  Monte  Carlo  Methods,  Lecture  Notes,  University  of  Bristol, 
Department  of  Mathematics,  Bristol,  England,  2007. 

9.  Oden,  J.T.  and  Demkowicz,  L.,  Applied  Functional  Analysis,  2nd  Ed.,  Computational 
Mechanics  and  Applied  Mathematics,  CRC  Press,  Boca  Raton,  Florida,  1996. 

10.  Weiner,  N.,  “The  homogeneous  chaos”,  American  Journal  of  Mathematics,  Vol.  60,  No.  4, 
1938,  pp.  897-936. 

11.  Xiu,  D.  and  Kamiadakis,  G.E.,  “The  Weiner- Askey  polynomial  chaos  for  stochastic 
differential  equations”,  SIAM  Journal  of  Scientific  Computing,  Vol.  24,  No.  2,  2002,  pp.  619- 
644. 

12.  Sandu,  A.,  Sandu,  C.  and  Ahmadian,  M.,  “Modeling  multibody  systems  with  uncertainties. 
Part  I:  Theoretical  and  computational  aspects”,  Multibody  System  Dynamics,  Vol.  15,  2006,  pp. 
373-395. 

13.  Soize,  C.  and  Ghanem,  R.,  “Physical  systems  with  random  uncertainties:  chaos 
representations  with  arbitrary  probability  measure”,  SIAM  Journal  of  Scientific  Computing,  Vol. 
26,  No.  2,  2004,  pp.  395-410. 


27 

Distribution  A 


14.  Hosder,  S.  and  Walters,  R.W.,  “Non-Intrusive  Polynomial  Chaos  Methods  for  Uncertainty 
Quantification  in  Fluid  Dynamics”,  AIAA  2007-129,  48th  AIAA  Aerospace  Sciences  Meeting,  4- 
7  January  2010,  Orlando,  Florida. 

15.  Hosder,  S.,  Walters,  R.W.  and  Balch,  M.,  “Efficient  Sampling  for  Non-Intrusive  Polynomial 
Chaos  Applications  with  Multiple  Uncertain  Input  Variables”,  AIAA  2007-1939,  48th  AIAA/ 
ASME/ASCE/AHS/ASC  Structures,  Structural  Dynamics  and  Materials  Conference,  23-26  April 
2007,  Honolulu,  Hawaii. 

16.  Abramowitz,  M.  and  Stegun,  I.A.,  Eds.,  Handbook  of  Mathematical  Tables  with  Formulas, 
Graphs  and  Mathematical  Tables,  Dover  Publications,  Mineola,  New  York,  1965. 

17.  Knio,  O.M.  and  Le  Maitre,  O.P.,  “Uncertainty  propagation  in  CFD  using  polynomial  chaos 
decomposition”,  Fluid  Dynamics  Research,  Vol.  38,  No.  9,  2006,  pp.  616-640. 

18.  Eldred,  M.S.  and  Burkardt,  J.,  “Comparison  of  Non-Intrusive  Polynomial  Chaos  and 
Stochastic  Collocation  Methods  for  Uncertainty  Quantification”,  AIAA  2009-976,  47th 
Aerospace  Sciences  Meeting,  5-8  January  2009,  Orlando,  Florida. 

19.  Burden,  R.L.,  Faires,  J.D.  and  Reynolds,  A.C.,  Numerical  Analysis,  2nd  Ed.,  Prindle,  Weber 
&  Schmidt,  Boston,  Massachusetts,  1981. 

20.  Walters,  R.W.,  “Towards  Stochastic  Fluid  Mechanics  Via  Polynomial  Chaos”,  AIAA  2003- 
0413,  41st  Aerospace  Sciences  Meeting,  6-9  January  2003,  Reno,  Nevada. 

21.  Hosder,  S.,  Walters,  R.W.  and  Perez,  R.,  “A  Non-Intrusive  Polynomial  Chaos  Method  for 
Uncertainty  Propagation  in  CFD  Simulations”,  AIAA  2006-891,  44th  AIAA  Aerospace  Sciences 
Meeting,  9-12  January  2006,  Reno,  Nevada. 

22.  Giunta,  A. A.,  Wojkiewicz,  Jr.,  S.F.  and  Eldred,  M.S.,  “Overview  of  Modern  Design  of 
Experiments  Methods  for  Computational  Simulations”,  AIAA  2003-649,  41st  Aerospace 
Sciences  Meeting,  6-9  January  2003,  Reno,  Nevada. 

23.  Wong,  T.-T.,  Luk,  W.-S.  and  Heng,  P.-A.,  “Sampling  with  Hammersley  and  Halton  Points”, 
Journal  of  Graphics  Tools,  Vol.  2,  No.  2,  1997,  pp.  9-24. 

24.  Cheng,  H.  and  Sandu,  A.,  “Collocation  Least-Squares  Polynomial  Chaos  Method”,  High 
Perfonnance  Computing  Symposium,  April  2010,  Orlando,  Florida. 

25.  Cheng,  H.,  Uncertainty  Quantification  and  Uncertainty  Reduction  Techniques  for  Large 
Scale  Simulations ,  Doctoral  Dissertation,  Virginia  Polytechnic  Institute  and  State  University, 

July  15,  2009,  Blacksburg,  Virginia. 

26.  L’Ecuyer,  P.,  “Uniform  random  number  generation”,  Annals  of  Operations  Research,  Vol. 

53,  1994,  pp.  77-120. 


28 

Distribution  A 


27.  Anderson,  J.D.,  Modern  Compressible  Flow:  With  Historical  Perspective,  2nd  Ed.,  McGraw- 
Hill,  New  York,  New  York,  1990. 

28.  Debusschere,  B.J.,  Najm,  H.N.,  Pebay,  P.P.,  Knio,  O.M.,  Ghanem,  R.G.  and  Le  Maitre,  O.P., 
“Numerical  challenges  in  the  use  of  polynomial  chaos  representations  for  stochastic  processes”, 
SIAM  Journal  of  Statistical  Computing,  Vol.  26,  No.  2,  2004,  pp.  698-719. 

29.  Field,  R.V.,  Jr.  and  Grigoriu,  M.,  “On  the  accuracy  of  the  polynomial  chaos  approximation”, 
Probabilistic  Engineering  Mechanics,  Vol.  19,  2004,  pp.  65-80. 


29 

Distribution  A 


APPENDIX  A 


DIRECT  STATISTICAL  ANALYSIS  FOR  TEST  PROBLEM  1 

This  problem  is  designed  to  be  a  simple  test  case  for  the  Legendre  non-intrusive 
polynomial  chaos.  Possessing  a  constant  integration  measure,  the  calculation  of  statistical 
parameters  is  greatly  simplified  requiring  only  the  evaluation  of  basic  integral  averages.  This  test 
case  involves  a  random  process  /  operating  on  a  single  random  variable  c  e  [-1,1]  .This  random 
process  is  written  as 


/(#)= i-jlf+f+f*) 


For  this  random  process,  the  probability  measure  is 

<Z)dZ  =X-dt; 

The  average,  or  expected  value,  for  /  is  computed  as  follows. 


Ue +£'+£*) 


By  evaluating  integrals,  we  obtain 


The  mean  for  this  random  process  is  therefore 

244 

E(f)= -  «  0.7746. 

315 


*-3 


j/'e  3  e5 


*  +  —  +  ■ 


(A-l) 


(A-2) 


(A-3) 


(A-4) 


(A- 5) 


For  this  problem,  the  computation  of  the  variance  is  similar.  Recall  that  the  definition  of  the 
variance  is 


var  (/)=  a2  =E(f2)-(E(f))2 


The  first  tenn  in  (A-6) 


£(/s)Aj 


Algebra  can  be  employed  to  simplify  the  integrand  obtaining 


(A-6) 


(A-7) 
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E(f )  =  A  1 1"1  -  A  (6?  +  5^4  +  4^ 6  -  3£8  -  2£10  -  £12  )1  cU;  (A-8) 

2 9 


This  integral  may  be  evaluated  as  follows. 


e3  e5  e7  e9  ell  el3  1 

6— +  5  — +  4- — 3— — 2— - (A-9) 

3  5  7  9  11  13 

-i 


By  simplifying,  this  moment  is 

E(f2)  =  1-^3  +  ----  —  -—]  *  0.66896  (A- 10) 

9^  7  3  11  13 

Substitution  of  (A-7)  and  (A- 5)  into  (A-6),  the  variance  can  be  computed  as 

var  (/)  =  a2  =  0.66896  -  (0.7746)2  =  0.068949  (A-l  1) 

Accordingly,  the  standard  deviation  is 

a  =  0.26258  (A- 12) 


E<f2)-L2i2)-Ts 
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APPENDIX  B 


DIRECT  STATISTICAL  ANALYSIS  FOR  TEST  PROBLEM  2 

This  problem  is  also  a  simple  test  case  for  the  Legendre  chaos,  but  it  is  designed  to 
incorporate  a  transcendental  function  into  the  random  process.  Again,  the  random  process 
operates  on  a  single  random  variable  d,  e  [-1,1] ,  and  it  is  defined  as 

m  =  l-£2)[sin(10£>  +  3]  (B-l) 

The  mean  of  this  random  process  may  be  computed  as  follows. 

E(f)  =  {7(1  -  £2)[sin(10£>  +  3]±-d£  (B-2) 

J,4  2 

E(f)  =  ^j  (1  - #2 )[sin(l 0|)  +  3] cU;  (B-3) 

8  1 

As  one  can  see,  the  probability  measure  remains  the  same  as  in  the  previous  case.  Equation  (B-3) 
is  evaluated  as  shown  below. 


£(/>A 


1  11  1 

J sin(10^)^  +  3 -  fg2  sin(10 £)<*£  -  3j  %2dd, 


(B-4) 


The  third  integral  in  (B-4)  can  be  evaluated  by  performing  a  detailed  integration  by  parts  as 


sin(10£)^  =  0  (B-5) 


Intuitively,  this  integral  may  be  easily  computed  by  realizing  that  the  integrand  is  anti¬ 
symmetric,  the  product  of  odd  and  even  functions.  For  this  reason,  it  must  equal  zero.  Hence, 


£(/)A 


j  sin(10^)^  +  3  -  3  j  fd% 


By  integrating,  we  have  that 


Eif)-\ 


cosqog) 

1  r\  - 


(B-6) 


(B-7) 
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£(/)  4 


cos(10)  |  2  |  cos(10)  |  2 


10 


10 


1 

2 


(B-8) 


By  using  similar  means,  the  variance  can  be  calculated.  First,  the  second  moment  of/  is 
determined,  i.e., 


E(f2)  =  j-^(l  -  2£2  +  ^4)(sin2  (10£>  +  6sin(10#)  +  9)\cU; 
16  1 


(B-8) 


By  expanding  the  integrand,  we  have  that 
1  rr 

E{f2)  =  —  [[sin2  (10/)  +  6sin(10/)  +  9  -  2£2  sin2  (10/)  - 12/  sin(lO^) 
32  J 


■18/+/  sin2  (1 0/  +  6/  sin(l  0/  +  9/ 


(B-9) 


With  a  significant  amount  of  effort,  the  integrals  in  (B-9)  can  be  evaluated  by  using  multiple 
stages  of  integration  by  parts  and  the  half-angle  formulas.  The  result  is 


9  1 

E(f2)  =  — 
32 


152  _  1657  .  ,  3 

- H - sin(20)  + - cos(20) 

15  40000  20000 


=  0.31785 


(B-10) 


Recalling  the  variance  formula 

var  (f)  =  E(f2)-(E(f))2  (B-ll) 

The  variance  is  computed,  through  (2-8),  (2-10)  and  (2-11)  as 

var(/)  =  a2  =  0.31785  -  (1/2)2  =  0.0678  (B-12) 

As  a  result,  the  standard  deviation  becomes 

cr  =  0.26048  (B-13) 
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APPENDIX  C 


DIRECT  STATISTICAL  ANALYSIS  FOR  TEST  PROBLEM  3 

This  problem  provides  a  two-dimensional  validation  test  case  for  the  Non-Intrusive 
Polynomial  Chaos  method  presented  earlier  in  this  report.  It  involves  a  non-linear  function  g  of 
two  random  variables,  i.e., 


g(x1?  x2)  =  ln(l  +  x12)sin(5x2)  (C-l) 

where  xi  and  a'2  are  uniform  random  variables  with  common  mean  2.0  and  a  common  coefficient 
of  variation  (CoV)  of  20%.  [C-l]  For  convenient  reference,  the  CoV  is  defined  as  the  ratio  of  the 
standard  deviation  to  the  mean,  i.e., 


CoV  =  —  (C-2) 

M 

These  uniform  random  variables  are  defined  on  the  interval  [a,  b\.  As  a  result,  x\  and  xi_  have  a 
common  probability  density  function  given  by 

f(x)  =  -r—  (c-3) 

b-a 

with  the  mean  (or  expected  value)  [C-2]  of 


M  = 


b  +  a 
2 


(C-4) 


The  variance  for  this  distribution  is  expressed  as 

(b-a)2 

12 

With  the  use  of  these  formulas,  it  may  be  shown  that 

a  =  2.0  -0.4-V3 
b  =  2.0  +  0.4a/3 


(C-5) 


(C-6) 


C.l  Expected  Value  of  g(xi,  jc2) 

Based  upon  this  information,  a  first  goal  is  the  compute  the  expected  value  of  g;  defined 
in  general  as 
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(C-6) 


oo  oo 

E(g(x \,x2))  =  J  { g(xl,x2)f(x1,x2)dx1dx2 

—00  —00 


where  /  is  the  joint  probability  distribution  function  of  x\  and  X2.[C-2]  These  two  random 
variables  are  independent,  so  the  joint  probability  function  is  calculated  as  the  product  of  their 
individual  probability  density  functions.  Hence, 


f(x  i,x2)  =  - — ^ 


1 


{b-aY  1.92 


(C-7) 


Equations  (g)  and  (h)  may  be  used  to  show  that 

b  b  /  i  S 

E\g(xx,  x2)]  -  [  fln(l  +  x2)  sin(5x2)  -  dxxdx2 

J  J  l 1.92 J 


Independence  and  the  product  fonn  of  g  can  be  used  to  rewrite  (C-8)  as 

1  b  b 

E[g(x  i,x2)]  =  — Jlna+xf  )dx]  | sin(5 x2)dx2 

a  a 

The  integrals  in  (C-9)  can  be  evaluated  independently,  i.e., 

b 

|  In  (1  +  x^)dxl  =  (xj  In  (1  +  xf)  -2xj  +  2  tan_1(Xj)| 


(C-8) 


(C-9) 


(C-10) 


b  j 

J sin (5x2  )dx2  =  —  [cos(5 a)  -  cos(5  b )] 


(C-ll) 


With  the  use  of  equations  (C-10)  and  (C-ll),  the  expected  value  of  g(x15x2)  is  calculated  as 


/j  =  E[g(xvx2 )]  =  0.079169 


(C-12) 


C.2  Variance  of  g(X|,  X2) 

Calculating  the  variance  of  g(x,  ,x2 )  is  substantially  more  complicated  than  the  expected 
value  computation.  The  variance  a2  of  a  random  variable  is  defined  as  follows.  [C-2] 

ct2  =E[{X-Mf-\  =  E[X2l[-M2  (C-13) 

To  obtain  the  variance  of  g,  it  is  necessary  to  obtain  the  moment  E[g2(xx,  x2)] .  That  is, 
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1  b  b 

E\g2(x  „x2)]  =  J  J  (ln(l  +  xf))  sin2 (5 x2)dxldx2 

^  *  a  a 


Fortunately,  the  multiple  integral  above  is  separable  in  that 


1  b  ^  b 

E[g2  (X, ,  x2 )]  =  — - -y  [  (ill  (1  +  x2  )fdxt  f  sin2  (5  x2 )  dx2 

(b~a)  Ja  { 


The  integral  over  the  random  variable  x2  is  easily  evaluated  as 

b  j  "  j 

[sin2(5  x2)dx2  =—  b-a  + — (sin(lOa)-sin(lOh)) 
J  ~  2  10 


(C-14) 


(C-15) 


(C-16) 


The  remaining  integral  is  much  more  difficult  to  evaluate.  In  fact,  it  cannot  be  evaluated  in  terms 
of  elementary  functions.  With  a  substantial  amount  of  work,  it  can  be  shown  that 


(C-17) 


J  (in (1  +  xf  ))<ix1  =  ^ (in (sec2  6*))"  tan  9  -  4  (ln(sec2  9)  -  2)(tan 9-9)  +  89 ln(cos  9) 

a 

eb 

-8j  k^cos#)^# 


where 


a  =  t&nOa  =>  0a=tan*(a) 
b  =  tan  0h  =>  0h  =  tan-1  (b) 


(C-18) 


The  final  integral  in  equation  (C-17)  cannot  be  evaluated  in  terms  of  elementary  functions.  It 
requires  the  use  of  numerical  quadrature  or  application  of  the  Clausen  function  CL2.  Appendix  E 
describes  the  Clausen  function,  but  the  core  result  needed  is  that 

e  , 

Jln(cos0)<70  =  —CL2(k -20)- 0\n2  ,  0  <9<-  (C-19) 

o  ^  ^ 

In  equation  (C-17),  the  corresponding  integral  may  be  written  as 

b  b  a 

J  In  (cos  9)d9  =  \  In  (cos  9)d9-\  In  (cos  9)  d9  (C-20) 

a  0  0 

Equations  (C-19)  and  (C-20)  may  be  combined  to  obtain  the  result 
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(C-21) 


h  1 

\\n(cosO)dO  =  -[CL2(x -26h)-CU{x -26a)\-(6h-6a)\n2 

ea  ^ 

Even  with  good  numerical  techniques  and  polynomial  series  expansions,  Clausen  functions  are 
challenging  to  evaluate.  For  this  reason,  the  definite  integral  in  equation  (t)  has  been  evaluated 
by  a  variety  of  techniques.  The  results  are  shown  in  Table  C-l.  For  the  numerical 
approximations,  the  composite  trapezoidal  and  Simpson’s  rules  are  applied.  Clausen  functions  1 
and  2  are  described  in  Appendix  E. 


Table  C-18.  Definite  Integral  of  In  (cos  6) 


Quadrature  Method 

Value 

Trapezoidal  Rule 

-0.2211678 

Simpson’s  Rule 

-0.2211675 

Clausen  Formula  1 

-0.2206295 

Clausen  Formula  2 

-0.2211681 

With  the  use  of  equations  (C-l 3),  (C-l 6)  and  (C-l 7)  with  the  data  in  Table  1,  variance 
approximations  can  be  computed  for  g(x i,  xi).  These  results  are  presented  in  Table  C-2. 


Table  C-2.  Variance  Approximations  for  Random  Function  g 


In  (cos  6)  Quadrature  Method 

Variance 

Standard  Deviation 

Trapezoidal  Rule 

1.263684 

1.124137 

Simpson’s  Rule 

1.263683 

1.124136 

Clausen  Formula  1 

1.262185 

1.123470 

Clausen  Fonnula  2 

1.263685 

1.124137 
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APPENDIX  D 


DIRECT  STATISTICAL  ANALYSIS  FOR  TEST  PROBLEM  4 

This  test  problem  involves  a  model  random  process  that  operates  on  four  uniform  random 
variables  x;  =xj{%j)  ,  j  =  1,...,4  .  Specifically,  the  process  is  written 

f(xl,x2,x3,x4)  =  exp[l  .5(x,  +x2  +x3  +  x4)]  (D- 

1) 

These  variables,  in  turn,  rely  upon  the  4  standard  uniform  variables  e [-1,1] .For  the 

random  variables,  x . ,  j  =  l,...,4,  the  mean  value  is  0.4  and  the  coefficient  of  variation  is  0.4.  [D- 
1]  From  this  information,  it  can  be  shown  that  the  domain  of  x  .  is 

x j  e  [a, b]  =  [0.4  -  0.08V12  ,0.4  +  0.08^12]  (D-2) 

Since  (D-l)  operates  on  four  independent  unifonn  random  variables,  the  joint  probability  density 
function  is  written  as 


w(x  l5x2,x3,x4)  = 


(b-a)4 


(D-3) 


The  expected  value  of  the  random  process  is  then  computed  as 

b  b  b  b  j 

£(/)=  1  [exp[l.5(x,  +  x-,  +  x3  +  x4 )] - -dx1  dx-i  dx3  dx4 


The  recurring  integral  appearing  in  (D-5)  is  easily  calculated  as 


J  exp(1.5  x)dx  =  £X^ 


exp(1.5x) 
5 


=  —  [exp(l  .5b)  -  exp(l  ,5n)] 


with  a  and  b  given  by  (D-2).  As  a  result,  the  expected  value  becomes 

2 


E(f)  = 


3 (b  -  a) 


(exp(l  .5b)  -  exp(l  .5n)) 


=  12.3609 


(D-4) 


j  b  b  b  b 

E(f)  - - - 1 cxp(l  .5x, )dxt  J exp(1.5x2)Jx2  J exp(1.5x3)Jx3  J exp(1.5x4)Jx4  (D-5) 

{b  ci)  a  a  a  a 


(D-6) 


(D-7) 


The  second  moment  of  this  random  process  may  be  calculated  in  the  same  way  after  noting  that 
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f\xvx2,x3,x4)  =  [exp(1.5(x1  +  x2  +  x3  +x4))]  ?=  exp[3(Xj  +x2  +  x3  +x4)]  (D-8) 


The  associated  integral  for  E(f2)  has  the  same  form  as  (D-5);  thus,  with 


b  ^ 

J  exp(3x)c£c  =  —  [exp(36)  -  exp(3o)] 


we  have  that 


E(f 2) 


1 

3  (b-a) 


(exp(3Z))  -  exp(3n)) 


(D-9) 


(D-10) 


Using  the  definitions  of  variance  and  standard  deviation,  the  standard  deviation  is  calculated  as 

cr  =  6.1556  (D-l  1) 
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APPENDIX  E 


CLAUSEN’S  FUNCTION 

In  association  with  the  polylogarithm,  polygamma  and  Riemann  Zeta  functions,  the 
Clausen  function  has  applications  in  physics,  such  as  in  the  field  of  electrodynamics.  [E-l]  As  is 
indicated  by  this  report,  it  also  has  uses  in  probability  modeling.  Clausen’s  function  is  defined  as 
follows. 


CL2  (0)  =  -  J  In 


e 

(  x\ 

hi 

2  sin 

X 

J 

0 

l2j 

dx 


(E-l) 


A  theoretical  series  expansion  also  exists  for  this  function,  i.e, 

Chi(e)=^pm  (e-2) 

i  k 

This  infinite  series  is  slowly  convergent  and  generally  not  used  for  calculations.  Instead,  series 
based  upon  Chebyshev  polynomials  are  often  recommended  for  generating  the  Clausen 
functions. [E-2, E-3]  Although  these  methods  are  accurate,  they  are  computationally  complicated 
and  still  involve  the  use  of  infinite  series.  These  difficulties  motivate  the  use  of  numerical 
quadrature  techniques  and  approximate  functional  forms. 

The  Clausen  integral  (E-l)  can  be  evaluated  by  the  use  of  basic  numerical  integration 
techniques  such  as  the  composite  forms  of  the  Trapezoidal  and  Simpson’s  rules. [E-4]  To  justify 
the  use  of  quadrature  on  (E-l),  it  is  necessary  to  show  that  the  integrand  is  finite  at  the  lower 
limit  of  integration.  At  first  glance,  the  sine  function  is  zero  at  zero,  so  the  natural  logarithm  at 
this  point  would  seem  to  tend  to  minus  infinity.  This  potential  difficulty  requires  more  in-depth 
analysis.  By  applying  the  asymptotic  approximation  for  the  sine  function  of  small  arguments,  we 
have  that 


sin 


x 

2 


In 


sin 


In' 


ln(x) 


(E-3) 


The  absolute  value  bars  have  been  eliminated  from  (E-3)  since  x  is  confined  to  the  domain  [0,  ji]. 
With  this  approximation,  we  can  estimate  the  Clausen  integral  as 


2  sin 


v 

dx  «  j  ln(x) dx  «  (x ln(x)  - . 


(E-4) 


If  this  result  is  finite  in  the  limit  as  x  — »  0 ,  then  the  Clausen  integral  is  finite  and  can  be 
approximated  by  quadrature.  This  limit  is  evaluated  as  follows. 
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Figure  E-l.  Approximate  values  of  the  Clausen  function  versus  an  angular  argument 


lim  xln(x) 

x->0 


lim 

jc-aO 


ln(x) 

1/ 

/  x 


lim(-x)  =  0 

x->0 


(E-5) 


L’Hopital’s  rule  has  been  employed  in  obtaining  this  result,  so  the  Clausen  integral  is  finite  and 
may  be  evaluated  by  use  of  the  composite  Trapezoidal  and  Simpson  rules.  Alternatively,  formula 
approximations  may  be  used  to  evaluate  the  Clausen  function. 

For  comparison,  there  are  two  approximate  formulas  for  the  Clausen  function. [E-5]  The 
simpler  formula  is  given  as 

CL2(0)  =  -01n (%)+  %4(k2  ~  92)-  {y2  -  21n2)sin<9  +  (in 2  -  1^6)sin2^  (E-6) 


where  0  <  9  <  rc .  The  second  formula,  regarded  to  be  of  higher  accuracy,  may  be  written  as 

CL2 (9)  =  -9  In sin(%)+  0/2%%O  ^  X120  ~  ^  +  W~ )  +  (2  ln 2  "  %)sin  9 

-  (89^28  -  In  2)sin  29  +  In  2  -  44^72)sin  3 9  (E-7) 

_  (425^/228g  -  ]/2  ln  2)sin  49  +  (% ln  2  -  1039/^75oo)sin  59 


where  0  <9  <n .  A  relative  error  of  0.63%  is  quoted  for  (E-6)  while  a  relative  error  of  0.003% 
is  similarly  quoted  for  (E-7).[E-3]  These  error  estimates  must  be  qualified  because  the  theory 
restricts  the  argument  9  to  a  rational  multiple  of  tt  where  a  rational  argument  is  the  ratio  of  two 
integers.  The  meaning  in  this  case  is  that  (E-6)  and  (E-7)  are  not  valid,  in  the  strictest  sense,  for 
irrational  multiples  of  n.  It  is  more  difficult  to  quantify  accuracy  in  this  case.  Still,  these  formulas 
seem  suitable  enough  for  performing  comparison  calculations.  A  plot  of  the  Clausen  function 
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evaluated  by  trapezoidal  rule  (TRAPZ),  Simpson’s  rule  (SIMPS),  formula  (E-6)  (APPRX1)  and 
formula  (E-7)  (APPRX2)  are  shown  in  Figure  (E-l). 

In  Appendix  C,  the  computation  of  variance  requires  evaluating  the  integral  of  the 
logarithmic  cosine,  i.e., 


|  ln(cosx)  dx 
0 

This  integral  can  be  evaluated  as  follows.  We  begin  with  the  Clausen  integral,  i.e., 

e 

CL2(6>)  = -Jln2sin|  —  ||  dx 


(E-8) 


6 

(v) 

hi 

2  sin 

X 

j 

0 

Iv 

(E-9) 


By  using  the  transformation  of  variables  x  =  2 y  along  with  the  trigonometric  identity 


sin 

fy  + 

=  2sin 

(y) 

cos 

(y) 

V2  2  j 

V2  j 

By  substituting  (E-10)  into  (E-9)  and  by  applying  a  property  of  the  logarithm, 


dy 


0/2 

en 

J  In 

2  sin 

i  y 

dy- 2  f  In 

2  cos 

i  y 

0 

UJ 

J 

0 

uJ 

2  cos 


fy) 

UJ 


dy 


Equation  (E-9)  easily  demonstrates  that 

en 

-CL,(60  =  -CL,(6>/2)-  f  In 
2  '  J 

With  the  change  of  variables  y  -  2x ,  we  have  that 

0/4  i  i 

|ln|2cosx|d'x  =  — CL2(6>)--CL2(0/2) 
o  4  2 

By  replacing  6  with  46* ,  we  obtain 

^  11 

J  ln|2cosx|  dx  =  — CL2(46>)  — CL2(20) 


(E-10) 


(E-l  1) 


(E-l  2) 


(E-l  3) 


(E-l  4) 


The  Clausen  identity 
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CL2(26>)  =  2CL2(<9)  -2CL2(;t  -20) 


(E-15) 


may  be  repeatedly  applied  to  obtain 

e  2 

J ln|2cosx|Jx  =  —CLfn  -20) 
0  2 


Hence, 


e 

J  ln|cosjc|  dx 

0 


^-CL2O-20)-01n2 


(E-16) 


(E-17) 


The  notation  CL2  denotes  a  Clausen  function  of  the  second  order.  Although  this  second 
order  function  is  more  commonly  used,  higher  order  Clausen  function,  e.g.,  CL4,  CL6,  exist. [E-3] 
These  functions  are  evaluated  from  series  of  orthogonal  polynomials.  This  numerical  procedure 
is  more  complicated  than  either  of  the  quadrature  rules  mentioned  earlier  or  the  formulas  above. 
The  methods  used  to  compute  the  CL2  create  curves  that  map  closely  together.  Hence,  we  have 
confidence  in  applying  quadrature  fonnulas  for  evaluating  Clausen’s  functions. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS  AND  ACRONYMS 


CFD . Computational  Fluid  Dynamics 

E . Error 

Hn . Hennite  polynomial 

LESLIE3D . Large  Eddy  Simulation  with  Linear  Eddy  Modeling  in  3  Dimensions 

Ln . Legendre  polynomial 

PC . Polynomial  chaos 

w . Probability  density  function 

a . Stochastic  or  random  process 

a,j . Expansion  coefficients 

<f„ . Random  variable 

^ . Vector  of  random  variables 

vFn . Orthogonal  polynomial 

9 . Random  event 

Q . Event  space 
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